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Abstract 

We  develop  a  deterministic  mathematical  model  to  describe  reactivation  of  latent  virus  by  chem¬ 
ical  inducers.  This  model  is  applied  to  the  reactivation  of  latent  KSHV  in  BCBL-1  cell  cultures 
with  butyrate  as  the  inducing  agent.  Parameters  for  the  model  are  first  estimated  from  known 
properties  of  the  exponentially  growing,  uninduced  cell  cultures.  The  model  is  then  extended  to 
describe  chemically  induced  KSHV  reactivation  in  latently  infected  BCBL-1  cells.  Additional 
parameters  that  are  necessary  to  describe  induction  are  determined  from  fits  to  experimental 
data  from  the  literature.  Our  model  provides  good  agreement  with  two  independent  sets  of 
experimental  data. 
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1  Introduction 


Many  viral  pathogens  establish  latency  and  are  dormant.  The  presence  of  inducers  leads  these  pathogens  to 
reactivate  and  replicate,  aiding  their  transmission  and  contributing  to  disease  development.  In  addition,  there 
is  increasing  evidence  in  the  literature  for  the  importance  of  polymicrobial  infections  in  which  microorganisms 
interact  in  a  synergistic  fashion,  impacting  both  pathogenesis  and  maintenance  of  health.  Among  these,  virus- 
bacteria  interactions  have  been  described,  including  reactivation  of  latent  virus  by  metabolic  end  products 
of  anaerobic  bacteria.  A  shift  in  the  balance  of  the  flora  often  controlled  by  the  intact  immune  system  may 
reflect  significant  morbidity  particularly  in  the  immune  suppressed  host.  The  relationships  between  viral 
pathogens  and  their  inducing  agents  have  not  previously  been  described  mathematically.  Therefore,  to  begin 
to  understand  the  relationship  between  pathogens  and  their  inducing  agents,  particulary  in  a  polymicrobial 
environment,  we  have  developed  a  mathematical  model  that  describes  the  reactivation  of  latent  herpes  virus 
by  metabolic  end  products  of  anaerobic  bacteria. 

Herpes  viruses  are  double-stranded  DNA  viruses.  Currently,  there  are  eight  known  herpes  viruses  that 
infect  humans.  After  primary  infection,  the  virus  remains  latent  in  specific  types  of  host  cells  that  may  be 
different  from  the  types  of  cells  targeted  for  primary  infection.  Latent  virus  persists  in  the  cell  nucleus  as 
episomal  DNA  until  it  is  reactivated,  beginning  a  program  of  lytic  replication  and  lysis,  and  leading  to  a  new 
(sometimes  asymptomatic)  round  of  infection  and  latency.  Upon  the  reactivation,  the  viral  DNA  begins  a 
lytic  program  characterized  by  a  temporal  cascade  of  gene  expression  culminating  in  excretion  of  free  virus 
and  cell  lysis.  The  lytic  program  is  typically  described  as  three  phases  of  gene  expression:  Immediate  Early, 
Early,  and  Late.  During  Immediate  Early  and  Early  phases  infected  cells  produce  viral  proteins  that  are 
necessary  for  viral  DNA  synthesis,  which  occurs  at  the  end  of  the  Early  phase.  During  the  Late  phase  of  the 
reproductive  cycle  the  host  cell  is  directed  to  make  the  structural  proteins  necessary  for  viral  packaging. 

The  exact  mechanisms  by  which  latent  virus  becomes  reactivated  and  begins  lytic  replication  are  not  entirely 
known.  However,  it  has  been  established  that  inducing  agents  such  as  Tetradecanoyl  Phorbol  Acetate 
(TPA),  sodium  butyrate,  and  other  short  chain  fatty  acids  (SFAs)  can  induce  lytic  replication  of  Kaposi’s 
Sarcoma-associated  Herpes  virus  (KSHV)  and  Epstein-Barr  virus  (EBV)  [14,  37,  38].  In  addition,  recent 
experiments  have  shown  that  the  spent  media  from  gram  negative  bacteria  cultures,  such  as  P.  Gingivalis 
and  P.  Intermedia,  which  contains  short  chain  fatty  acids  (e.g.,  iso- valeric,  n-butyric  acid,  and  propionic 
acid),  can  also  induce  latent  KSHV  to  begin  lytic  replication  [29]. 

A  mechanism  for  bacterial  reactivation  of  latently  infected  cells  has  strong  health  implications  for  the  oral 
environment  as  well  as  the  gut  and  GI  tract,  where  there  may  be  large  numbers  of  gram  negative  bacteria  in 
the  presence  of  latently-infected  cells.  Reactivation  of  latent  herpes  viruses  are  a  major  health  concern  for 
immune-compromised  individuals,  such  as  those  with  AIDS.  Understanding  the  role  of  anaerobic  bacteria  in 
reactivation  of  latent  herpes  viruses  may  have  important  health  consequences  if  there  is  similar  reactivation 
of  latent  episomal  HIV  by  bacterial  metabolic  end  products. 

To  our  knowledge  there  has  been  no  mathematical  modeling  treatment  of  viral  reactivation  at  the  cellular 
level.  Much  of  the  mathematical  modeling  of  herpes  viruses  has  focused,  understandably,  on  modeling  at  the 
epidemiological  level  (e.g.,  [7,  18]).  Recently,  Wang  et  al.,  used  HHV-6  infection  as  a  stimulus  for  studying 
cellular  changes  in  the  T  cell  immune  system  under  pathological  conditions  [34] .  Their  study  included  data 
from  the  literature,  clinical  data,  and  cell  culture  data.  Their  model  agreed  well  with  data,  but  focused  on 
the  viral  load  and  T  cell  response.  Clearly  there  is  a  need  and  opportunity  to  understand  viral  latency  and 
reactivation,  especially  since  this  is  a  characteristic  feature  of  herpes  virus  infection. 

In  this  manuscript  we  report  on  a  deterministic  mathematical  model  that  we  have  developed  to  describe 
reactivation  of  latent  virus  by  chemical  inducers.  In  particular,  we  apply  this  model  to  the  reactivation  of 
latent  KSHV  in  BCBL-1  cell  cultures  with  butyrate  as  the  inducing  agent.  KSHV,  also  known  as  Human 
Herpesvirus-8  (HHV-8),  is  a  gamma  herpes  virus  that  is  responsible  for  Kaposi’s  Sarcoma  tumor  development 
and  other  lymphoproliferative  disorders  such  as  Castleman’s  disease  and  primary  effusion  lymphoma.  KSHV 
latently  infects  epithelial  and  lymphoid  cells  that  are  present  in  the  oral  environment  and  reactivation  of 
latent  virus  in  B  cells  may  play  a  role  in  the  pathogenesis  of  Kaposi’s  sarcoma  [25].  BCBL-1  cells  are  an 
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immortalized  cell  line  derived  from  body-cavity-based  lymphomas  that  are  latently  infected  with  multiple 
copies  of  KSHV,  but  not  Epstein-Barr  virus  (EBV)  [28],  providing  a  convenient  vehicle  for  studies  of  viral 
reactivation.  We  choose  butyrate  as  the  inducing  agent  since,  as  an  SFA,  it  is  similar  to  the  metabolic 
end  products  of  gram  negative  bacteria,  it  is  easily  quantified,  and  its  activity  is  well  documented  in  the 
experimental  literature. 

We  first  estimate  parameters  for  our  model  from  known  properties  of  the  exponentially  growing,  uninduced 
cell  cultures.  We  then  extend  the  model  to  describe  chemically  induced  KSHV  reactivation  in  latently  infected 
BCBL-1  cells.  Additional  parameters  that  describe  induction  are  determined  from  fits  to  experimental  data 
available  in  the  literature.  Our  model  exhibits  good  agreement  with  two  independent  sets  of  experimental 
data. 

The  model  presented  here  establishes  a  general  framework  for  modeling  the  effect  of  other  inducing  agents 
that  act  through  histone-deacetylase  (HDAC)  inhibition,  including  other  SFA’s  produced  by  the  metabolic 
processes  of  gram  negative  bacteria.  As  such,  it  may  also  be  applied  to  other  latent  virus  systems  that  are 
induced  to  replicate  via  HDAC  inhibition,  such  as  EBV,  HIV,  and  HCMV. 


2  Modeling  Compartments 


We  describe  the  dynamics  of  the  host  cells  and  viral  DNA  copies  using  a  set  of  ordinary  differential  equations 
(ODEs).  A  schematic  diagram  of  the  ODE  model  compartments  is  shown  in  Figure  1  and  the  model 
compartments  are  summarized  in  Table  1.  Latent  L  and  lytic  R  copies  of  viral  DNA  reside  in  the  nuclei 
of  host  cells.  As  mentioned  above,  there  are  multiple  copies  of  episomal  viral  DNA  in  each  BCBL-1  host 
cell.  We  make  the  following  “all  or  nothing”  simplifying  assumption:  within  a  given  host  cell,  nuclear  viral 
DNA  copies  are  either  all  latent  or  all  in  a  lytic  replication  program.  Therefore,  there  are  two  types  of  host 
cells,  host  cells  Hr  with  lytic  virus  only  or  host  cells  Hr  with  latent  virus  only.  In  future  work,  we  may 
superimpose  a  probability  distribution  on  the  parameters  to  better  approximate  mixed  conditions  where  a 
host  cell  may  contain  both  latent  and  lytic  virus  in  varying  levels.  Such  a  modeling  technique  was  successfully 
used  in  cellular  level  HIV  models  to  account  for  variable  length  (with  uncertainty)  pathways  in  [2] .  In  models 
of  this  type  the  state  variables  are  the  expected  values  of  concentrations  (or  of  numbers  of  cells)  resulting  in 
delay  differential  equation  models  embodying  uncertainty.  We  do  not  pursue  this  level  of  refinement  in  the 
initial  model  developed  here. 

Host  cell  death  can  occur  as  a  result  of  natural  aging,  with  individual  rates  da  and  viral  lysis  rate  di 
[30],  or  the  inducing  agent  rates  Al(s)  and  d#(s).  Host  cells  that  die  are  added  to  a  nonviable  host  cell 
compartment  N.  During  the  Early  Phase  of  the  lytic  program,  productive  replication  of  viral  DNA  takes 
place  in  replication  compartments  [35] .  The  localization  of  the  replication  process  allows  for  exponential-like 
growth  of  viral  DNA  R ,  where  the  progeny  of  replicating  virus  become  replication  templates  themselves 
[16].  Capsid  assembly  occurs  in  the  nucleus  and  nucleocapsid  envelopment  occurs  at  the  nuclear  membrane. 
The  intracellular  viral  DNA  compartment  Vi  represents  viral  DNA  copies  that  are  no  longer  targets  for 
replication  and  are  available  for  encapsulation. 

After  envelopment,  the  virus  is  released  as  free  virions  Vp.  Experiments  by  Bechtel  et  a l.,  indicate  that  free 
KSHV  virions  produced  by  BCBL-1  cells  fail  to  infect  many  cultured  lymphoid  cell  lines  [5].  Therefore,  we 
assume  that  free  virions  that  are  produced  are  not  capable  of  reinfecting  the  host  cells.  We  can  also  arrive 
at  this  assumption  by  considering  that,  within  the  cultured  cell  line,  there  are  always  a  small  number  of 
host  cells  with  spontaneously  reactivated  lytic  virus  that  are  producing  free  virions.  If  the  free  virus  were 
able  to  reinfect  the  BCBL-1  cells  then  the  average  number  of  viral  copies  per  host  cell  would  not  be  a  stable 
quantity.  Since  there  is  no  indication  that  the  average  number  of  viral  DNA  copies  in  the  cell  line  is  changing 
with  time,  we  assume  that  free  virions  are  not  reinfecting  the  host  cells. 
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Figure  1:  Schematic  of  the  modeling  compartments  associated  with  latent  virus  host  cells  Hr,  lytic  virus  host 
cells  Hr,  and  nonviable  host  cells  N.  Latent  virus  L  reactivates  to  become  lytic  virus  R,  either  spontaneously 
or  in  response  to  an  inducing  agent  s.  Lytic  virus  R  undergoes  exponential  growth  until  it  passes  to  the 
intracellular  viral  compartment  Vj  where  it  is  available  for  packaging  and  is  released  as  free  virions  Vf- 


Compartment 

Symbol 

Units 

Host  cells  (latent  virus  only) 

hl 

number  of  cells 

Host  cells  (lytic  virus  only) 

Hr 

number  of  cells 

Nonviable  cells 

N 

number  of  cells 

Latent  virus 

L 

DNA  copies 

Lytic  Virus 

R 

DNA  copies 

Intracellular  virus 

Vr 

DNA  copies 

Free  virus 

Vf 

number  of  virions 

Table  1:  ODE  Model  Compartments. 


3  Mathematical  Model  for  the  Uninduced  Case 


Host  cell  dynamics.  Latently  infected  host  cells  Hr  are  growing  at  a  rate  PlHr  and  have  a  natural  death 
rate  cIlHl,  where  /3l  >  d-L-  As  part  of  the  lytic  cycle,  herpes  viruses  block  the  cell  cycle  in  Go/Gi  and  block 
cell-initiated  apoptosis  [17,  20,  36].  Therefore,  host  cells  with  lytic  virus  Hr  are  expected  to  have  much 
smaller  rates  of  cell  replication  /3rHr  and  natural  death  (IrHr  than  host  cells  with  latent  virus  (/3r  <C  Pl 
and  d,R  <C  dL).  Although  lytically  infected  host  cells  Hr  are  expected  to  have  a  reduced  death  rate  due  to 
viral  inhibition  of  apoptosis,  there  is  an  additional  mechanism  for  cell  death  due  to  the  production  of  virus 
and  the  resulting  cell  lysis.  We  model  this  viral-induced  death  rate  as  a  function  of  the  average  amount 
of  intracellular  virus  that  accumulates  cIi(Vi)Hr,  where  Vr  =  Vi/Hr.  Nonviable  cells  are  in  a  process  of 
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disintegration  into  smaller  fragments  and  leave  the  N  compartment  with  a  rate  pN. 


Experimental  observations  of  uninduced  cell  cultures  indicate  that  some  fraction  of  host  cells  will,  through 
spontaneous  reactivation,  have  lytic  virus  [8,  23,  27].  Therefore,  we  include  terms  for  spontaneous  reactivation 
of  latently  infected  cells  with  a  rate  constant  ao-  When  herpes  viruses  establish  latency  in  newly  infected 
cells  they  undergo  a  small  lytic  replication  program  using  proteins  imported  in  the  viral  capsid.  Rather  than 
proceeding  to  encapsulation  and  lysis,  the  replicated  copies  of  viral  DNA  then  enter  a  latent  state  [16].  This 
initial  lytic  infection  that  leads  to  latency  is  different  from  the  lytic  replication  that  occurs  when  latent  virus 
in  cultured  cells  is  induced  into  a  lytic  program  and  it  is  not  clear  whether  virus  that  is  reactivated  from 
latency  can  return  to  a  latent  state  in  these  cultured  cells.  Therefore,  we  allow  for  the  possibility  that  lytic 
virus  returns  to  a  latent  state  with  a  rate  pR.  Later  we  will  present  an  argument  that  this  term  must  be 
zero  in  a  cell  culture  with  a  stable  average  viral  copy  number.  With  these  considerations  in  mind,  we  model 
the  host  cell  dynamics  in  the  uninduced  case  by 


dHL 

dt 


(7l  -  a0)HL  +  pHR 


dHR 

dt 


otoHL  +  (7 R  ~  P~  di(Vi ))  H R 


dN 

dt 


dRHL  +  (dR  +  d;(Vi ))  Hr  —  pN , 


(3.1) 


where  7 l  =  Pl  -  dL  and  jR  =  j3R  -  dR. 


Viral  dynamics.  There  are  multiple  copies  of  viral  DNA  in  each  latently  infected  host  cell  nucleus,  with 
an  average  L  per  cell.  Latent  virus  is  maintained  as  circular  episomal  DNA  in  the  host  cell  nucleus  where  it 
is  tethered  to  the  host  cell  DNA  and  is  copied  each  time  the  host  cell  DNA  is  replicated.  Therefore  we  can 
approximate  the  L  compartment  growth  at  a  rate  which  is  proportional  to  the  host  cell  growth  rate  and  the 
average  number  of  latent  virus  copies  per  host  cell,  L/3lHr .  We  allow  for  the  possibility  that  lytic  virus  R 
is  also  copied  during  host  cell  replication,  although  as  stated  above,  we  expect  the  host  cells  with  lytic  virus 
to  be  replicating  at  a  much  lower  rate.  This  gives  a  similar  growth  rate  of  R/3rHr  for  lytic  virus. 


Using  the  same  reasoning,  we  assume  that  each  dying  host  cell  destroys  a  number  of  latent  (lytic)  virus 
equal  to  the  average  number  of  latent  (lytic)  virus  per  host  cell,  as  represented  by  the  loss  term  LdRHR 
(RdRHR  and  Rdi(Vi)HR).  We  assume  that  the  intracellular  viral  DNA  copies  are  not  replicated  with  the 
host  cell,  since  they  are  not  tethered  to  the  host  DNA,  but  are  destroyed  when  the  host  cell  dies  at  rates  of 
Vidj(Vi)HR  and  VidRHR.  The  reactivation  and  deactivation  of  latent  and  lytic  virus  follows  in  a  manner 
similar  to  the  terms  above.  Therefore,  the  uninduced  model  for  the  viral  DNA  dynamics  is  given  by  the  set 
of  equations 

dL  -  -  - 

,,  =  7  lLHl  —  a0LHL  +  pRHR 

dt 

(1f.  =  (K-q)R  +  a0LHL+  (■yR-dI(yI))HRR-  pHRR 


dU7 

dt 


qR  —  pVi  —  dRHRVi  —  di(Vi)HRVi 


(3.2) 


dVp 

dt 


pVi, 


where  kR  is  the  replication  rate  for  lytic  viral  DNA,  qR  is  the  rate  at  which  lytic  virus  moves  to  the 
intracellular  compartment,  and  pVj  is  the  rate  at  which  the  intracellular  DNA  is  packaged  and  excreted  as 
free  virions  Vf- 
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Uninduced  cell  and  viral  dynamics. 

full  uninduced  model  as 

dHL 

dt 

dHR 

dt 

dN 

dt 

dL 

dt 

dR 

dt 


Using  L  =  L/Hl,  R  =  R/Hr  and  VT 
(7 l  -  a0)HL  +  pHR 

&0HL  +  (7 R  ~  di  (V»  -  p)  Hr 

dhRlh  +  [dR  +  di(Vi ))  Hr  —  pN 

(7 l  ~  ol0)L  +  pR 

(k  -  q  +  7_r  -  di(Vi)  -  p)  R+a0L 


Vi/Hr  ,  we  can  write  the 


(3.3) 


dV7 

dt 


qR  -  (p  +  dR  +  d/(U/))  Vi. 


The  solution  for  Vp  is  easily  obtained  from  the  above  solutions  and  can  be  written  as 


VF(t)  =  VF0  +  [  pVi(u)du. 
Jt0 


(3.4) 


4  Mathematical  Model  for  the  Induced  Case 


We  next  modify  the  uninduced  model  (3.3)  to  include  the  actions  of  an  inducing  agent  s.  The  main  affect  of 
the  inducing  agent  is  to  increase  the  rate  at  which  latent  virus  becomes  reactivated  (a(s)  >  an)-  In  addition, 
inducing  agents  such  as  butyrate  and  valproate  may  also  cause  host  cell  death  through  activation  of  host  cell 
genes  [6,  22,  37].  We  represent  the  induced  cell  death  rates  by  SR(s)HR  and  Sl(s)Hl ,  where  5p(s)  and  SR(s) 
are  functions  of  the  butyrate  concentration  s.  With  these  additional  terms  the  equations  for  the  induced 
case  become 


and 


dHL 

dt 

dHR 

dt 

dN 

dt 

dL 

dt 

dR 

dt 

dVi 

dt 


(7 l  -  a(s)  -  SL(s ))  Hl  +  pHR 

(7 R  ~  dR(s)  —  di(Vi ))  Hr  +  a(s)HL  —  pHR 

( dL  +  Sl(s))Hl  +  ( dR  +  SR(s)  +  di(Vi ))  HR 

(7 l  -  <*(s)  -  h(s))  L  +  pR 

(k  -  q  +  7fl  -  SR(s)  -  di(Vi))  R  +  a(s)L  -  pR 

qR  -  (p  +  dR  +  diiVi)  +  8r(s))  Vi 


VF{t)  =  VF0  +  f  pVi(u)du. 
Jt0 


(4.5) 


(4.6) 
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5  Parameter  Estimation 


Some  of  the  parameters  in  the  uninduced  model  (3.3)  can  be  estimated  from  physiological  considerations. 
The  parameter  q  is  the  rate  constant  for  movement  of  lytically  replicating  virus  from  the  R  compartment  to 
the  intracellular  compartment  Vj.  We  estimate  that  the  rate  at  which  lytic  virus  leaves  the  R  compartment 
is  inversely  proportional  to  the  time  that  the  lytic  virus  spends  in  the  lytic  program  of  gene  expression  and 
replication.  Therefore,  we  set  q  =  1  /T,  where  T  is  the  approximate  time  it  takes  to  complete  the  lytic 
program  (approximately  48  hours). 

Since  the  lytic  virus  causes  cell  cycle  arrest  and  inhibits  apoptosis,  we  make  the  simplifying  assumptions  that, 
for  lytically  infected  host  cells  HR,  the  growth  rate  constant  and  natural  death  rate  constant  both  vanish, 
i.e.,  7_r  =  dR  =  0.  Although  lytically  infected  host  cells  HR  are  expected  to  have  a  reduced  rate  of  death 
due  to  viral  inhibition  of  apoptosis,  there  is  an  additional  mechanism  for  cell  death  due  to  the  production 
of  virus  and  the  resulting  cell  lysis.  We  have  modeled  the  rate  for  this  viral-induced  death  as  a  function  of 
the  average  amount  of  intracellular  virus  that  accumulates  which  is  denoted  by  dj(Vr),  where  Vj  =  Vi/Hr. 
For  the  sake  of  simplicity  and  in  the  interest  of  obtaining  parameter  estimates,  we  assume  di  to  be  a  linear 
function  of  the  average  number  of  intracellular  copies  per  lytically  infected  host  cell,  i.e.,  dj{Vj)  =  cVj. 

Implications  of  the  p  term  on  cell  culture  stability.  We  expect  that  the  average  number  of  latent 
copies  of  viral  DNA  per  latently  infected  host  cell  is  a  constant  that  is  characteristic  of  the  particular  cell 
line;  otherwise  the  cell  properties  would  be  different  from  one  experiment  to  the  next.  Therefore,  L/ HR  =  L 
is  a  constant,  which  we  designate  with  the  variable  n ,  and  dL/dt  =  n{dHR/dt).  Using  these  conditions  and 
the  first  and  fourth  equations  in  (3.3),  we  find  that  p{R  —  nHR)  =  0.  This  condition  can  be  met  in  two 
ways,  either  p  =  0  or  R  =  nHR.  The  former  condition  requires  that  there  is  no  reversion  of  lytically  infected 
host  cells  back  to  latency.  The  latter  condition  implies  that  dR/dt  =  n(dHR/dt)  and,  from  the  second  and 
fifth  equations  in  (3.3),  this  condition  would  require  that  k  —  q  =  0.  Since  kR  is  the  rate  of  replication  of 
lytic  virus  and  qR  is  the  rate  at  which  the  lytic  virus  is  moving  to  the  intracellular  virus  compartment,  the 
condition  n  —  q  =  0  implies  that  there  can  be  no  net  accumulation  of  replicating  virus.  If  k  —  q  =  0  there 
will  be  no  exponential  growth  of  viral  DNA  copies.  Such  a  condition  is  contrary  to  current  knowledge  of  the 
lytic  cycle  for  KSHV  [16]. 

In  summary,  if  the  condition  L  =  n  is  to  be  maintained  in  the  uninduced  BCBL-1  cell  cultures  then  either 
there  can  be  no  reversion  of  virus  back  to  latency  (p  =  0)  or  the  rate  at  which  the  lytic  viral  DNA  is 
being  copied  must  be  equal  to  the  rate  at  which  viral  DNA  copies  are  being  added  to  the  Vj  compartment 
(k  —  q  =  0).  The  latter  condition  implies  that  copies  of  replicated  viral  DNA  can  not  be  templates  for 
replication,  themselves,  which  contradicts  the  current  understanding  of  the  lytic  replication.  Therefore,  we 
set  p  =  0  in  this  model.  It  should  be  noted  that  these  conditions  only  apply  in  the  case  of  stable,  latently 
infected  cell  lines  such  as  BCBL-1  cells  and  do  not  apply  to  acutely  infected  cells  where  lytic  virus  is  known 
to  revert  to  latency. 

Properties  of  cell  lines.  Other  parameters  in  the  uninduced  model  (3.3)  can  be  estimated  from  known 
properties  of  the  uninduced  cell  cultures.  As  we  have  noted,  BCBL-1  cells  are  an  immortalized  cell  line 
derived  from  body-cavity-based  lymphomas  that  are  latently  infected  with  multiple  copies  of  KSHV  [28]. 
Exponentially  growing  cells  are  maintained  in  a  medium  consisting  of  growth  nutrients  and  antibiotics  to 
prevent  bacterial  contamination.  Growing  cells  are  periodically  split  to  prevent  contact  growth  inhibition 
and  to  provide  fresh  nutrients.  We  expect  that,  under  constant  growth  and  maintenance  conditions,  certain 
average  properties  of  the  cell  line  will  be  unchanging  in  time.  Some  of  the  parameters  for  the  uninduced 
model  can  be  estimated  from  these  constants. 

Table  2  summarizes  these  constants,  along  with  ranges  of  values  reported  in  the  experimental  literature.  One 
such  constant  is  the  fraction  of  host  cells  with  spontaneously  reactivated  virus  as.  Experimental  observations 
of  this  quantity  vary  from  1-8%.  Differences  in  observed  values  for  as  may  be  due  to  the  different  methods 
used  to  detect  lytic  virus  as  well  as  differences  in  the  growth  and  maintenance  conditions  of  the  cell  lines. 
Another  constant  is  the  fraction  of  nonviable  cells  Nr  in  the  cell  culture,  with  values  reported  from  8-20%. 
Variability  in  the  observed  fraction  of  nonviable  cells  Nr  is,  also,  most  likely  a  result  of  differences  in  cell 
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line  growth  and  maintenance  conditions  as  well  as  differences  in  measurement  techniques  ( e.g .,  dye  exclusion 
versus  dye  inclusion  identification  techniques  and  haemocytometer  versus  cell  sorting  counting  techniques). 
BCBL-1  cells  have  a  doubling  time  Dp  of  approximately  48  to  85  hours,  which  will  also  be  dependent  upon 
growth  and  maintenance  conditions.  The  model  parameter  n,  which  represents  the  average  number  of  latent 
virus  copies  per  uninduced  host  cell  is  a  constant  of  the  cell  line,  independent  of  growth  and  maintenance 
conditions.  This  quantity  is  not  measured  directly,  but  rather  a  related  quantity  n^,  representing  the  average 
number  of  copies  of  viral  DNA  per  cell  (or  average  cell-associated  viral  DNA,  i.e. ,  L  +  R  +  Vj),  is  measured. 
Values  of  nx  ranging  from  50  to  70,  dependent  upon  growth  and  maintenance  conditions,  have  been  reported 
in  the  literature. 

In  Table  2  the  uninduced  cell  line  constants  are  specified  in  terms  of  the  ODE  model  compartments.  The 
subscript  t  — >  oo  indicates  that  we  expect  these  quantities  to  be  constants  in  an  equilibrated  average  sense 
under  conditions  of  constant  growth  and  maintenance.  In  addition  to  the  above  observed  cell  line  constants, 
we  also  expect  that  the  average  number  of  lytic  Ra  and  intracellular  Via  DNA  copies  per  host  cell  to  be 
constant  though  we  have  no  experimental  measurements  for  these  constants. 


Constant  Description 


Symbol  Value  Units  Model  Formulation 


Fraction  of  lytic  host  cells  as 

[13,  26,  33] 

Fraction  of  nonviable  host  cells  Nr 

[37,  38] 


Host  cell  doubling  time 
[24,  33] 


D„ 


Average  copies  of  viral  DNA  per  cell  nx 
[23,  27] 

Average  number  of  lytic  DNA  per  cell  Ra 

Average  number  of  intracellular  Via 

DNA  per  cell 


0.01-0.08 


0.08-0.2 


50-70 


Hr 


Hl  +  Hr +  N 


Nr  = 


N 


HL  +  Hr +  N 


48-85  hr  2  = 


(Hl  +  Hr  +  N)t=Dp 
(Hl  +  Hr  +  N)t=o 


tit  = 


L  +  R  + Vi 
Hl  +  Hr  +  N 

R 

Hr  /  t—>oo 


Hr 


Table  2:  Observed  cell  line  constants. 
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From  the  constants  and  relationships  in  Table  2,  the  uninduced  equations  (3.3),  and  the  above  estimates  for 
q ,  (in,  7#,  and  di(Vj)  =  cVi,  we  obtain  the  following  parameter  values  in  terms  of  the  unknown  parameters 


7  L,  H, 


dr.  = 


'  1L 


ln(2  )/Dp  +  pNr 
1  —  as  —  Nr 

~^—{Nr  -  1)  +  (1  -  a,  -  Nr) 

asVIADpK  ’  asVIAK  ’ 


«0  —  7  L 


ln(2) 

Dn 


P  = 


Ra  (1  -as-Nr)(  ln(2)\ 


TVia  a. 

7 l  ~  ln(2 )/Dp 


7  l 


Dr 


-) 


D  i  (Ra  —  NrRA  —  ut  +  asVjA)  + 
asRA  J  1 


n  = 


nr  —  as(RA  +  Via) 
1  —  as  —  Nr 


(5.7) 


6  Numerical  Results 


In  this  section,  we  compare  our  reactivation  model  to  two  sets  of  experimental  data  from  the  literature  that 
describe  reactivation  of  latent  virus  in  BCBL-1  cells.  In  particular,  we  compare  our  model  to  longitudi¬ 
nal  cell  viability  data  following  chemical  induction.  In  one  case,  Zoeteweij,  et  al.,  [38]  use  n-butyric  acid 
(CH3CH2CH2COOH)  as  the  inducing  agent,  while  in  the  other  case,  Yu,  et  al.,  [37]  use  Na  n-butyrate 
(C H3C H2C H2COO N a)  as  the  inducing  agent.  Despite  the  differences  in  molecular  formulas,  we  expect  the 
two  inducers  to  behave  roughly  the  same  in  solution  after  dissociation  of  the  H+  and  Na+  ions.  In  particular, 
for  cell  viability  data,  Sakurazawa,  et  al.,  have  shown  that  n-butyric  acid  and  Na  n-butyrate  are  cytotoxic 
at  approximately  the  same  concentrations  [31]. 


6.1  Parameter  values 

In  Section  5  we  estimated  values  for  some  parameters  from  physiological  considerations  and  showed  that 
other  model  parameters  could  be  estimated  from  a  few  free  parameters  and  experimental  observations  of 
cell-line  constants  (5.7).  The  fraction  of  nonviable  cells  Nr  (in  the  uninduced  case)  is  the  only  value  that 
we  can  determine  directly  from  the  data  of  Zoeteweij,  et  al. ,  and  Yu,  et  al.,  obtaining  approximate  values  of 
8%  and  17%,  respectively.  For  the  other  constants  we  choose  values  from  within  the  ranges  reported  in  the 
literature  (Table  2):  as  =  0.017,  Dp  =  85  hr.  In  addition,  we  choose  rir  =  70  and  63  for  modeling  the  data 
of  Zoeteweij,  et  al.,  and  Yu,  et  al..,  respectively.  In  this  way  we  have  n  =  74  for  both  sets  of  data.  We  choose 
the  following  values  for  the  free  parameters  and  unknown  constants  :  7 l  =  8.4  x  10-3  hr-1,  /x  =  2.1  x  10~4 
hr-1,  Ra  =  89,  and  Via  =  89. 

In  Table  3  we  tabulate  the  parameter  values  used  in  simulations  corresponding  to  the  above  choices  of 
experimental  constants  and  free  parameters. 


Parameter 

Symbol 

Zoeteweij,  et  al.,  data 

Yu,  et  al.,  data 

Units 

Net  growth  rate  constants  for  host  cells 

7  L 

8.4  x  1CT3 

8.4  x  10-3 

hr-1 

7  R 

0 

0 

hr-1 

Natural  death  rate  constant  for  host  cells 

di 

6.49  x  10"4 

1.70  x  10"4 

hr"1 

0 

0 

hr-1 

Spontaneous  reactivation  rate  constant 
for  latently  infected  host  cells 

a0 

2.45  x  10"4 

2.45  x  10"4 

hr-1 

Spontaneous  deactivation  rate  constant 
for  lytically  infected  host  cells 

P 

0 

0 

hr"1 

Rate  constant  for  cell  death  due  to  viral  lysis 

c 

5.48  x  10"5 

3.99  x  10"5 

hr"1 

Rate  constant  for  nonviable  cell  degradation 

p 

2.1  x  10"4 

2.1  x  10"4 

hr-1 

Rate  constant  for  synthesis  of  viral  DNA 

K 

2.30  x  10"2 

2.28  x  10"2 

hr"1 

Rate  constant  for  sequestration  of 
viral  DNA  for  encapsulation 

q 

2.08  x  10"2 

2.08  x  10"2 

hr"1 

Rate  constant  for  packaging  and 
secretion  of  virions 

p 

7.80  x  10"3 

9.13  x  10"3 

hr-1 

Average  number  of  copies  of 
viral  DNA  per  latent  host  cells 

n 

74 

74 

- 

Induced  reactivation  rate  constant 

OLc 

8.72  x  10_1 

1.85  x  10_1 

hr-1 

Induced  death  rate  constant 

Sc 

5.33  x  10"3 

6.91  x  10"3 

hr-1 

Induced  death  rate 

SL(s) 

0 

0 

hr"1 

Table  3:  Parameters  from  the  uninduced  model  (3.3)  are  calculated  from  (5.7)  with  constants  as  =  0.017, 
Via  =  89,  Ra  =  89,  Nr  =  0.08  or  0.17,  Dp  =  85  hr,  and  nr  =  70  or  63,  and  with  free  parameters 
=  8.4  x  10^3  hr-1  and  /r  =  2.1  x  10-4  hr-1.  Parameters  from  the  induced  model  (4.5)  are  obtained  from 
fits  to  experimental  data. 


6.2  Uninduced  Case 

We  first  present  results  of  a  simulation  for  the  uninduced  case  modeled  by  system  (3.3).  Parameter  values 
used  for  the  simulation  are  shown  in  Table  3.  The  initial  condition  for  all  compartments  is  zero,  except 
for  compartments  Hr  and  L ,  which  have  initial  conditions  of  1.0  x  106  and  7.4  x  10'.  Figure  2  depicts  the 
fraction  of  nonviable  cells  N/(Hr  +  Hu  +  N)  (dashed  lines)  and  the  fraction  of  spontaneously  reactivated 
cells  Hu/ (Hl  +  Hu  +  N)  from  simulations  using  equations  (3.3)  and  the  parameter  values  in  Table  3.  In 
Fig.  2  we  can  see  that,  by  1000  hours,  the  fraction  of  nonviable  cells  and  the  fraction  of  spontaneously 
reactivated  cells  asymptotically  reach  the  specified  equilibrium  values  of  as  =  0.017  and  Nr  =  0.08  for  the 
data  of  Zoeteweij,  et  al.,  and  as  =  0.017  and  Nr  =  0.17  for  the  data  of  Yu,  et  al.  By  equilibrium  we  mean 
that,  although  the  cell  culture  is  growing  exponentially,  certain  characteristic  properties  related  to  the  ratios 
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of  model  compartments  eventually  become  constants  (see  Table  2) . 


(a) 


(b) 


Figure  2:  Uninduced  simulations  using  equations  (3.3)  and  parameters  from  Table  3  for  the  data  of  a)  Zoeteweij,  et 
al,  and  b)  Yu,  et  al.  Solid  lines  plot  the  percentage  of  cells  that  are  spontaneously  reactivated  Hr/(Hl  +  Hr  +  N), 
while  dashed  lines  plot  the  percentage  of  cells  that  are  nonviable  cells  N/(Hl  +  Hr  +  N). 


In  Fig.  3  we  can  see  a  similar  equilibration  of  the  average  number  of  lytic  and  intracellular  DNA  copies 
per  lytically  infected  host  cell,  R/Hr  and  Vi/Hr ,  respectively.  In  Fig.  3  we  can  see  that,  by  1000  hours, 
the  quantities  R/Hr  an  Vi/Hr  have  reach  the  specified  equilibrium  values  of  Ra  =  89  and  Via  =  89, 
respectively  (see  Table  2). 


(a)  (b) 

Figure  3:  Uninduced  simulations  using  equations  (3.3)  and  parameters  from  Table  3  for  the  data  of  a)  Zoeteweij, 
et  al.,  and  b)  Yu,  et  al.  The  plots  show  the  average  number  of  lytic  viral  DNA  copies  R/ Hr  (dashed  line)  and  the 
average  number  of  intracellular  viral  DNA  copies  Vi/Hr  (solid  line)  per  lytically  infected  host  cell.  Both  quantities 
start  with  an  initial  condition  of  zero,  but  the  quantity  R/ Hr  grows  quite  rapidly,  giving  the  appearance  of  having 
a  nonzero  starting  value  for  these  plot  limits. 


The  equilibrated  simulations  for  the  uninduced  model  approximate  the  properties  of  the  uninduced  cell 
cultures  that  are  subsequently  used  in  induction  experiments  and  provide  initial  conditions  for  simulation  of 
the  induced  equations  (4.5). 
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6.3  Induced  Case 


Next  we  report  on  simulations  for  the  induced  equations  (4.5).  In  the  case  of  those  parameters  that  are 
common  to  both  the  uninduced  and  induced  models  we  use  the  same  values  as  in  the  previous  simulation 
(Table  3).  Initial  conditions  for  the  induced  model  are  obtained  from  the  simulations  of  the  uninduced 
model  at  t  =  2400  hr.  Not  only  is  the  inducing  agent  capable  of  reactivating  the  latent  virus,  but  it  can 
also  activate  host  cell  genes  that  may  lead  to  apoptosis.  We  make  a  simplifying  assumption  that  activation 
of  host  cell  genes  would  most  likely  occur  in  conjunction  with  reactivation  of  latent  virus.  In  other  words, 
we  assume  that  inducing  mechanisms  that  would  initiate  cell  apoptosis  would  also  lead  to  viral  reactivation 
and,  therefore,  assume  that  there  will  be  no  chemically  induced  cell  death  for  latently  infected  host  cells 
(i.e.,  5l  =  0). 

The  exact  functional  forms  of  the  rates  for  induced  lytic  cell  death  5r(s)  and  induced  reactivation  of  latent 
cells  cr(s)  are  not  known.  We  first  choose  simple  affine  functions  a(s)  =  ao  +  acs  and  <5r(s)  =  6cs  and  find 
values  for  the  function  parameters  by  fitting  longitudinal  experimental  data  on  BCBL-1  cell  viability  from 
Zoeteweij,  et  al.,  [38]  and  Yu,  et  al. ,  [37],  separately.  The  parameter  fitting  is  accomplished  by  forming  an 
ordinary  least  squares  inverse  problem  as  described  in  the  Appendix  and  then  estimating  the  parameters 
using  a  Nelder-Mead  algorithm.  Standard  errors  are  calculated,  the  details  of  which  are  also  given  in  the 
Appendix.  In  the  case  of  the  affine  functions,  estimated  values  for  the  parameters  are  insensitive  to  the 
initial  values  that  seed  the  optimization  algorithm. 

Figure  4  compares  cell  viability  predicted  from  simulations  with  data  from  both  experimental  groups,  using 
the  estimated  parameters  for  Sc  and  ac  obtained  by  the  ordinary  least  squares  estimation  techniques.  Es¬ 
timated  parameter  values  are  reported  in  Table  3  and  4.  Some  of  the  model  parameters  differ  between  the 
two  groups  because  of  differences  in  uninduced  cell  viability  for  the  two  groups  (92%  versus  83%).  From 
Fig.  4  it  can  be  seen  that  the  simulations  for  the  induced  model  qualitatively  match  the  behavior  of  the 
experimental  data. 


(a) 


(b) 


Figure  4:  Comparison  of  cell  viability  measurements  with  simulations  using  induced  equations  (4.5)  and  fitted 
parameters  for  linear  functions  a(s)  and  <5r(s):  a)  Zoeteweij,  et  al,  circles  0  mM,  squares  0.03  mM,  triangles  0.3 
mM,  and  diamonds  3  mM,  ac  =  0.872,  Sc  =  5.33  x  10_  !  and  b)  Yu,  et  al,  circles  0  mM,  triangles  0.3  mM,  stars  1.5 
mM,  and  diamonds  3  mM,  ac  =  0.185,  =  6.91  x  10_,i. 
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In  Fig.  5  we  plot  the  number  of  free  virions  as  a  function  of  time  for  different  butyrate  concentrations.  Yu, 
et  al. ,  [37]  observe  in  their  experiments  that  high  concentrations  of  butyrate  (1.5  and  3  mM)  greatly  increase 
lytic  activity,  but  also  significantly  increase  cell  death.  The  end  result  is  that,  even  after  5  days,  very  few 
free  virions  are  produced  because  of  massive  amounts  of  cell  death  before  the  end  of  the  lytic  program.  This 
is  contrasted  by  observations  at  smaller  concentrations  of  butyrate  (<  0.3  mM),  where  much  less  cell  death 
is  seen  and  there  is  significant  secretion  of  free  virions.  In  Fig.  5  it  can  be  seen  that  there  is  approximately 
a  three-fold  increase  in  free  virion  produced  at  0.3  mM  concentration  of  butyrate  as  compared  to  the  3  mM 
concentration. 


(a) 


(b) 


Figure  5:  The  results  of  simulations  of  free  virions  produced  using  induced  equations  (4.5)  and  parameters  for  linear 
functions  a(s)  and  5r(s)  fitted  to  experimental  data  from  a)  Zoeteweij,  et  al. ,  and  b)  Yu,  et  al. 


Table  4  summarizes  the  estimated  parameters,  standard  errors,  and  confidence  intervals  obtained  from  fitting 
the  induced  equations  (4.5)  with  the  parameter  functions  a(s)  =  «o  +  acs  and  Sr(s)  =  Scs  to  both  sets  of 
data.  Table  4  shows  that  the  estimated  parameter  values  for  both  groups  are  within  an  order  of  magnitude 
of  each  other.  Differences  between  the  parameter  values  may  reflect  differences  in  the  cell  growth  and 
maintenance  conditions  or  differences  in  experimental  measurement  techniques.  Even  in  the  uninduced  case, 
there  is  a  difference  in  the  cell  viability  for  both  groups,  with  Nr  =  0.08  for  the  data  of  Zoeteweij,  et  al.,  and 
Nr  =  0.17  for  the  data  of  Yu,  et  al.  In  addition,  Zoeteweij,  et  al,  measure  cell  viability  using  Dead  Red 
staining  and  flow  cytometry,  while  Yu,  et  al.,  measure  cell  viability  with  trypan  blue  staining  and  counting 
on  a  haemocytometer.  In  Table  4,  it  can  also  be  seen  that  the  standard  errors  for  the  reactivation  rate 
constants  ac  are  at  least  an  order  of  magnitude  less  than  the  parameter  values.  However,  the  standard 
errors  for  the  induced  death  rate  constants  Sc  are  the  same  order  of  magnitude  as  the  parameter  values, 
providing  us  with  less  confidence  in  values  obtained  for  6C. 
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Data 

Parameter 

Estimated  Value 

Standard  Error 

Confidence  Interval 

Zoeteweij,  et  al. 

<5C 

OLc 

5.33  x  1CT3 

8.72  x  1CT1 

3.91  x  10"3 

1.26  x  10~2 

[-3.11  x  10-3, 1.38  x  KT2] 
[8.45  x  10_1, 8.99  x  10"1] 

Yu,  et  al. 

Sc 

OLc 

6.91  x  10"3 

1.85  x  HT1 

3.70  x  10"3 

1.11  x  10"2 

[-1.03  x  10— 3 , 1.48  x  10"2] 
[1.61  x  10_1, 2.09  x  10"1] 

Table  4:  Estimated  parameter  values,  standard  errors,  and  confidence  intervals 


7  Discussion 


In  other  simulations,  we  used  different  functional  forms  for  a(s)  and  5r(s),  including  Michaelis-Menton  and 
sigmoid  functions,  but  we  found  that  the  fits  of  the  induced  equations  to  cell  viability  data  were  relatively 
insensitive  to  more  complicated  functional  forms  (data  not  shown)  and  that  reasonable  fits  to  cell  viability 
data  were  obtained  by  assuming  simple  linear  functions.  However,  with  more  data,  especially  with  data  for 
viral  DNA  compartments,  we  expect  to  be  able  to  determine  optimal  functional  forms  for  a(s)  and  <5r(s), 
for  example,  combinations  of  linear,  Michaelis-Menton,  and  sigmoid  functions.  Alternatively,  instead  of 
fixing  the  functional  form  of  a(s)  and  (5fj(s)  a  priori  in  parametric  form,  we  could  estimate  the  shape  of  the 
functional  form  itself  using  approximation  by  piece-wise  linear  splines  or  other  approximations  as  has  been 
successfully  done  in  other  problems  in,  for  example,  [1,  3]. 

Even  though  this  preliminary  model  yields  good  qualitative  agreement  with  cell  viability  data,  additional 
experimental  data  is  needed  to  compare  model  predictions  to  other  measurable  quantities  of  interest,  such 
as  cell-associated  DNA  (L  +  R  +  Vi)  and  free  virions  (VV).  Additional  data  would  also  help  to  determine 
the  free  parameters  7 l  and  /r,  as  well  as  the  unknown  constants  Ra  and  Via-  In  the  case  of  Ra  and  Via , 
parameter  sensitivity  tests  show  that  the  optimal  parameter  values  Sc  and  ac  are  relatively  insensitive  to 
variations  in  Ra  or  Via,  since  varying  Ra  or  Via  by  ±  5%  produced  3%  or  less  variation  in  the  optimized 
parameter  values. 

Instead  of  a  single  viral  compartment  R  to  quantify  copies  of  viral  DNA  in  the  lytic  program,  we  could 
modify  the  model  to  describe  Immediate  Early,  Early,  and  Late  gene  expression  (RNA),  represented  by 
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compartments  Ri,  R2,  and  R3  in  the  following  equations  for  an  induced  model 


and  Vp(t) 


dHL 

dt 

dHn 

dt 

dN 

dt 

dL 

dt 

dR\ 

dt 

dR2 

dt 

dR 3 
dt 

dV, 

dt 


(7 l  -  a(s)  -  SL(s ))  Hl  +  pHR 

(7 r  ~  $r(s)  -  d/(Vr)  -  p)HR  +  a(s)HL 

{cLl  +  Sl(s))Hl  +  (dR  +  SR(s)  +  di(Vi ))  HR  —  pN 

(7 l  -  a(s)  -  <5l(s))  L  +  p(Ri  +  R2  +  R3) 

(7 R  ~  Qi  -  £r(«)  -  ^/(Vz)  -  p)  f?i  +  a(s)L 
(k  +  7fl  -  92  -  <5fl(s)  -  dj(Vr)  -  p)  i?2  +  9i^i 
(lR-Q-  Sr(s)  -  d^Vi)  -  p)  R3  +  q2R2 
qR3  -  (p  +  dfl  +  d/(V»  +  6r(s))  Vr 


^fo  +  [  pVi{u)du. 
Jto 


(7.8) 


The  three  parameters  q\,  q2l  and  q3  represent  the  rate  at  which  viral  DNA  moves  from  one  stage  of  the  lytic 
program  to  the  next.  These  parameters  can  be  estimated  as  l/Tf,  I/T2,  and  I/T3,  respectively,  where  T\, 
T2,  and  T3  are  the  approximate  times  for  each  stage  of  the  lytic  program.  Corresponding  parameters  in  this 
proposed  model  and  (4.5)  would  not  necessarily  represent  the  same  quantities. 


By  having  model  compartments  that  quantify  RNA  production  or  promoter  activity  from  genes  representative 
of  each  stage  of  the  lytic  cycle,  we  can  hope  to  predict  viral  reactivation  in  more  detail  and  compare 
to  experimental  gene  expression  data.  For  example,  ORF50,  vIL6,  K8.1  could  be  representative  of  the 
Immediate  Early,  Early,  and  Late  stages,  respectively.  A  single  compartment  L  can  represent  latent  gene 
expression,  primarily  ORF73  expression. 


There  may  be  underlying  biological  delays,  due  to  the  ordered  cascade  of  gene  expression  that  makes  up  the 
lytic  program,  that  are  not  captured  with  the  model  (4.5).  A  model  such  as  (7.8)  in  which  we  rewrite  the 
single  R  compartment  as  three  compartments  R±,  R2,  and  R3  representing  the  Immediate  Early,  Early,  and 
Late  phases  of  the  lytic  program  might  be  expected  to  more  closely  capture  the  biological  delays  inherent  in 
the  lytic  program  of  system. 


8  Conclusion 


We  have  developed  a  preliminary  deterministic  mathematical  model  to  describe  reactivation  of  latent  virus 
by  chemical  inducers.  In  particular,  we  apply  this  model  to  the  reactivation  of  latent  KSHV  in  BCBL-1  cell 
cultures  with  butyrate  as  the  inducing  agent.  We  first  estimate  parameters  for  our  uninduced  model  from 
physiological  considerations  and  known  properties  of  these  exponentially  growing,  uninduced  cell  cultures. 
We  then  extend  the  model  to  describe  chemically  induced  KSHV  reactivation  in  latently  infected  BCBL-1 
cells.  Additional  parameters  that  describe  induction  are  determined  from  fits  to  experimental  data  available 
in  the  literature.  Our  model  provides  good  agreement  with  two  independent  sets  of  experimental  data. 
While  this  preliminary  model  yields  good  qualitative  agreement  with  cell  viability  data  for  KSHV  induced 
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by  butyrate,  it  also  strongly  suggests  the  need  for  further  experiments  designed  explicitly  to  support  model 
development  and  validation  in  providing  not  only  more  but  also  additional  types  of  longitudinal  data. 

This  model  exploits  the  polymicrobial  nature  of  reactivation  by  focusing  on  chemical  inducers  that  utilize  the 
same  mechanisms  as  gram  negative  anaerobic  bacteria.  Among  viruses  there  can  be  interactions  that  result 
in  reactivation  as  well.  For  example,  recent  experiments  have  shown  that  KSHV  enhances  HIV  replication 
and  reactivation  [9,  11]  and  that  KSHV  ORF50  (lytic)  gene  products  increase  in  vitro  cell  susceptibility 
to  human  immunodeficiency  virus  type  1  infection  [10].  This  evidence  for  the  synergistic  interactions  of 
KSHV  and  HIV  highlights  the  importance  of  further  model  development  and  application  to  polymicrobial 
environments. 

Application  of  the  type  of  model  we  have  developed  here  to  other  latent  viruses  could  provide  information 
about  the  similarities  and  differences  among  latent  virus  systems  and  could  define  the  relationship  of  these 
organisms  to  their  inducers.  For  example,  like  KSHV,  many  latent  viruses  are  induced  to  replication  via 
HDAC  inhibition  and  are  responsive  to  agents  like  sodium  butyrate.  These  include  EBV,  HCMV,  HSV, 
HIV,  Adenovirus,  HPV,  and  HTLV1.  The  model  with  appropriate  variations  could  potentially  be  applied 
to  any  of  these  viruses.  Further,  the  use  of  this  and  future  models  with  other  inducers  may  also  provide 
extremely  valuable  clinical  information  about  induction  mechanisms. 


9  Appendix 


In  this  appendix  we  discuss  the  asymptotic  theory  used  to  compute  the  standard  errors  and  confidence 
intervals  in  Table  4  of  Section  6.3.  We  first  give  a  general  summary  of  the  theory. 

We  assume  N*  scalar  longitudinal/inducer  level  observations  (time/inducer  series  of  numbers  or  ratios  of 
numbers  of  cells  as  described  below)  are  represented  by  the  statistical  model 

Yj  =  fj(60)  +  ej,  j  =  1,2, . . .  N* ,  (9.9) 

where  fj(9 o)  is  the  model  for  the  observations  in  terms  of  the  state  variables  and  9q  £  Rm  is  a  “set”  of 
theoretical  “true”  parameter  values  (assumed  to  exist  in  a  standard  statistical  approach).  We  assume  for 
our  statistical  model  of  the  observation  or  measurement  process  (9.9)  that  the  errors  €j,  j  =  1,2, ... ,  N* , 
are  independent  identically  distributed  ( i.i.d .)  random  variables  with  mean  E[ej]  =  0  and  constant  variance 
var[ej]  =  <7q,  where  of  course  ctq  is  unknown  (standard  residual  plots  with  the  data  used  in  our  simulation 
suggested  this  assumption  of  constant  variance).  We  then  have  that  the  observations  Yj  are  i.i.d.  with  mean 
E[Yj ]  =  fj(9o)  and  variance  var[Yj\  =  erg. 

We  consider  estimation  of  parameters  using  an  ordinary  least  squares  (OLS)  approach.  Thus  we  seek  to  use 
data  {yj}  for  the  observation  process  {Yj}  with  the  model  to  seek  a  value  9  that  minimizes 


JV* 

j(9)  =  Y.\y,-fM\2-  (9-10) 

3  = 1 

Since  Yj  is  a  random  variable,  we  have  that  the  estimator  9ols  is  also  a  random  variable  with  a  distribution 
called  the  sampling  distribution.  Knowledge  of  this  sampling  distribution  provides  uncertainty  information 
(e.g.,  standard  errors)  for  the  numerical  values  of  9  obtained  using  a  specific  data  set  { yj}  (i.e.,  a  realization 
of  {Fj})  when  minimizing  J{9). 

Under  reasonable  assumptions  on  smoothness  and  regularity  (the  smoothness  requirements  for  model  so¬ 
lutions  are  readily  verified  using  continuous  dependence  results  for  ordinary  differential  equations  in  our 
example;  the  regularity  requirements  involve,  among  others,  conditions  on  how  the  observations  are  taken  as 
sample  size  increases,  i.e.,  N*  — ■>  oo),  the  standard  nonlinear  regression  approximation  theory  ([15],  [19],  [21], 
and  Chapter  12  of  [32])  for  asymptotic  (as  N*  — >  oo)  distributions  can  be  invoked.  This  theory  yields  that 
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the  sampling  distribution  9(Y)  for  the  estimate  9,  where  Y  =  {Yj}^=1,  is  approximately  a  m-multivariate 
Gaussian  with  mean  E[9(Y)\  and  covariance  matrix  cov[0(Y)\  =  S0  =  <Xq[xT  (9o)x{^o)]~1  ■  Here  y(0)  =  Fg{9) 
is  the  N*  x  m  sensitivity  matrix  with  elements 

XM  =  and  Fe{9)  =  (fie(9), . . . ,  fN.e(0))T. 

That  is,  for  N*  large,  the  sampling  distribution  approximately  satisfies 

Ools{Y)  ~  Afm(e0,a20[xT(e0)x(9o)}-1)  :=  Um(9o,Y0).  (9.11) 


The  elements  of  the  matrix  \  =  (Xjk)  can  be  estimated  using  the  forward  difference 

XlM  -~m T  ~ - h, - ■ 

where  h k  is  an  m- vector  with  nonzero  entry  in  only  the  kth  component,  or  using  sensitivity  equations  (see 
[4]  and  the  references  therein).  For  our  efforts  here  we  chose  the  sensitivity  equation  approach  as  explained 
below.  Since  9q,  <j$  are  not  known,  we  must  approximate  them  in  So  =  Co[xT (0o)x(@o)]~1  ■  For  this  we 
follow  standard  practice  and  use  the  approximation 

s0«s  (6)  =  &2ixT(e)x(e)}-1 

where  6  is  the  parameter  estimate  obtained,  and  the  approximation  cf2  to  er2  is  given  by 

CTo  «  v2  =  N*^m  Y  I  Vi  ~  fi$)  I2' 
j'=i 


Standard  errors  to  be  used  in  confidence  interval  calculations  are  thus  given  by  SE^(B)  =  yS kk{9),  k  = 
1, 2, . . . ,  m  (see  [12]). 


In  the  induced  case  example  of  Section  6.3,  we  consider  the  parametric  functional  forms  6r(s)  =  dcs  and 
a(s)  =  acs  +  «o-  If  we  let  x  =  (Hl,  Hr,N,  L,  R,Vi,Vf)T  and  denote  9  =  ( 6c,ac ),  then  the  differential 
equations  in  the  induced  case  can  be  written  in  a  general  form 


x  =  g(t,x,s,0) 

ir(0)  =  a;0, 


(9.12) 


where  g  :  K+  xl"  x  R+  x  — >  Mn  for  n*  =  7,  m  =  2,  and  Xq  =  ( Hl0 ,  Hr  0,  N0,  L0,  R0,Viq,  Vfo)T  ■  Since 

the  experimental  data  are  given  in  percentage  of  viable  cells,  we  define  the  outputs  of  the  model 


Vtotal{t,S,9)  -  N(t,  s,  9) 
Ytotal  (t,  S,  $) 


t,S>  0, 


where  Vtotai  =  Hr  +  Hr  +  N.  In  each  parameter  fit,  we  use  data  that  is  longitudinal  (taken  at  tk )  and  across 
several  levels  s,;  of  inducer.  This  is  indexed  by  Tj  =  (tk,  s-t)  for  k  =  1, . . . ,  K,  i  =  1, . . . ,  /,  and  observations 
Dj  for  the  model  values  fj{9)  =  f(tk,  Si,  9),j  =  1, . . . ,  N*  =  KI.  Then,  we  construct  the  OLS  estimator  by 
minimizing  the  cost  criterion  (9.10)  where  {yj}  denotes  the  experimental  data  (in  the  data  of  Section  6.3  we 
had  N*  =  15  or  16  resulting  from  K  =  4  and  7  =  4-  see  Figure  4).  For  the  optimization  in  9  we  used  the 
Nelder-Mead  algorithm. 


dF  * 

To  compute  the  covariance  matrix  S  we  need  the  sensitivity  matrix  Fg.  That  is,  y(0)  =  -^r— (0).  From  the 

o9 

dx 

outputs  defined  in  (9.12),  it  suffices  to  have  the  sensitivities  — .  To  compute  these  we  used  the  sensitivity 
equation  method  which  involves  solving  the  n*  x  m  matrix  variational  differential  equation 


d  /  dx  \  dg  dx  dg 
Jt\d9J  =  d^d9+  W 


(9.13) 
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where  the  matrix  coefficient  and  the  forcing  function  in  this  equation  are  evaluated  along  solutions  of  the 
system  equation  (9.12).  Note  that  this  variational  equation  can  be  solved  simultaneously  (see  [4]  for  details) 
with  the  system  equation  (9.12). 

Finally,  in  order  to  compute  the  confidence  intervals  (at  the  100(1  —  c)%  level)  for  the  estimated  parameters 
in  our  example,  we  define  the  confidence  level  parameters  associated  with  the  estimated  parameters  so  that 

P(0k  -  tc/2SEk(§ )  <ek<ek+  tc/2SEk0))  =  1  -  C,  (9.14) 

where  c  £  [0,1],  and  tc/2  £  R+.  For  a  given  c  value  (small,  say  c  =  .05  for  95%  confidence  intervals),  the 
critical  value  tc/2  is  computed  from  the  Student’s  t  distribution  tN  ~m  with  N*  —  m  degrees  of  freedom 
since  for  each  of  the  data  sets  available  to  us  we  have  N*  is  less  than  30.  The  value  of  tc/2  is  determined  by 
P{T  >  tc/2)  =  c/2  where  T  ~  tN*~m. 
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